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Abstract 



We study the dynamic critical behavior of the multi-grid Monte Carlo (MGMC) 
algorithm with piecewise-constant interpolation applied to the two-dimensional 0(4)- 
symmetric nonlinear a- model [= SU{2) principal chiral model], on lattices up to 
256 x 256. We find a dynamic critical exponent j^i = 0.60 ± 0.07 for the W- 
cycle and z int j^p. = 1.13 ± 0.11 for the V-cycle, compared to z int j^i = 2.0 ± 0.15 for 
the single-site heat-bath algorithm (subjective 68% confidence intervals). Thus, for this 
asymptotically free model, critical slowing-down is greatly reduced compared to local 
algorithms, but not completely eliminated. For a 256 x 256 lattice, W-cycle MGMC is 
about 35 times as efficient as a single-site heat-bath algorithm. 
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1 Introduction 



By now it is widely recognized ]I|, 0, that better simulation algorithms, with 

strongly reduced critical slowing-down, are needed for high-precision Monte Carlo studies 
of statistical-mechanical systems near critical points and of quantum field theories (such as 
QCD) near the continuum limit. One promising class of such algorithms is multi-grid Monte 
Carlo (MGMC) |5], || |7|]. In paper I of this series [f§, we explained the conceptual founda- 
tions of the MGMC method, and proved the absence of critical slowing-down for Gaussian 
(free-field) MGMC. In paper II [[?[], we applied MGMC to the two-dimensional ferromagnetic 
XY model, and analyzed the dynamic critical behavior in the high-temperature (vortex) 
and low-temperature (spin- wave) phases (see also ||). 

In the present paper we apply MGMC to the two-dimensional 0(A r )-symmetric nonlinear 
cr-model (also called iV-vector model) with N = 4, defined by 



H(a) 



-f3o{4) 



(xx>) 



where each spin a x is a unit vector in R 4 , (xx') runs over all nearest-neighbor pairs, and 
/?o(4) > is the inverse temperature in the 0(4) normalization. Since the sphere S 3 is 
isometric (as a Riemannian manifold) to the group SU(2) via the map 
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a + ia 3 



a 2 + ia 1 



1.2) 



a 



ia 



this model can also be thought of as the SU(2) principal chiral model 

H(g) = -Psu{i) Retr(4fe') , 

{xx 1 ) 



1.3) 



with Psu{2) = /3o(4)/2- Our principal aim is to make a careful study of the dynamic critical 
behavior of the MGMC algorithm for this asymptotically free model, and in particular to 
extract the dynamic critical exponent (s). We also take advantage of the high efficiency of 
the MGMC algorithm — about 35 times that of a single-site heat-bath algorithm, on a 
256 x 256 lattice — to make a detailed comparison of the static scaling behavior with the 
asymptotic-freedom predictions. 

Two-dimensional iV-vector models are of direct interest in condensed-matter physics, 
and they are important "toy models" in elementary-particle physics by virtue of their close 
similarity to four-dimensional gauge theories. They have therefore been extensively studied 
by rigorous |5|, |It], [IT], [T2|, [TJ| , exact-but-nonrigorous ]TJ], [T5|, [H], [T7|, |TB|, [Tj| , renormalization- 

group pg, p, n, n, m H 13 13, n, m 



high-temperature expansion 30, 31 



1/JV 

expansion [^, g§ |5j, f37|, [BJ and Monte Carlo 0, |4jJ, g |§ |§ |5J g7|, |§ 
methods.^ It is known rigorously that every infinite-volume Gibbs measure is invariant 
under global spin rotations fly]] and in particular has zero spontaneous magnetization ; it 



1 This is a very incomplete set of references, but a significant fraction of the relevant articles can be tracked 
down starting from these. Also, studies dealing solely with the XY case are cited in and therefore omitted 
here. 
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is believed, though not yet proven, that at each (3 there is a unique infinite-volume Gibbs 
measure .0 

Renormalization-group calculations in the low-temperature expansion (= weak-coupling 
perturbation theory) ^UJ ^TJ suggest that the iV-vector models with N > 2 are asymptotically 
free — i.e. that the only critical point is at ft = oo — and that the correlation length £ and 
susceptibility x behave as 



= Q ^-1/(^-2)^/9/(^-2) 



1 + ^ + ^ + 

ft ft 2 



1 + ^ + ^ + 

ft P 2 



'M) 

1.5) 



as P — > oo.f] Here the and are nonuniversal constants (depending on N) that can be 
computed in weak-coupling perturbation theory on the lattice at k + 2 loops; in particular, 
ai and b% have been computed for the nearest-neighbor action ( |1.1[ ) and two variant actions 
p6| , |29|.P| On the other hand, and C x are nonuniversal constants (depending on N) that 
have to be computed by some non-perturbative method. For the exponential correlation 
length (= inverse mass gap), the coefficient Q = 2~ 5 / 2 [(N - 2)/(2ne n / 2 )] 1 / ( - N - 2 \m/A m )- 1 
hasQ been computed exactly using the Bethe Ansatz |17], [Tj| [I - 



the result is 



m 



l/(7V-2) 



r i + 



JV-2 



1.6) 



sec 



54, 55, 56 



The asymptotic-freedom predictions (|1.4j ) / (|1.5|) have not yet been proven rigorously (but 
501, [H], Q for some progress), and indeed they have been questioned by one group 
However, they are consistent with Monte Carlo studies for N 
]V = 4 @H § and Af = 8 @. More precisely, for N 
P < 2.05 (corresponding to 10 < f < 300) 



45, 48 



3 it is found that for 1.5 ^ 
the quantity defined by ( |1.4|) varies by about 
20% — contrary to the prediction that it should be constant — and also differs by about 
20% from the predicted exact value ( |1.6| ). It appears, therefore, that quantitative asymptotic 
scaling has not yet set in at £ ~ 100, even though the qualitative behavior is correct. For 



2 In the XY case it has been proven that there is a unique translation-invariant infinite-volume Gibbs 
measure [[l2| p^| . 

3 In this paper the susceptibility is defined as x — Sx( cr o ■ c x ). This definition is employed also in 
{H |(| |H ||, [H fto|, Q. By contrast, the standard thermodynamic definition is Xthermo = Pxthis paper] 
the extra factor of arises because the inverse temperature (3 is introduced as an overall factor multiplying 
also the magnetic field. The thermodynamic definition is employed in ]2l| , 28 
must be borne in mind when comparing results in the literature! 

4 Reference pT 
equation (12a) of 



This difference in notation 



employs, as we do, the second-moment definition of the correlation length: compare 
§ with (p^)/(|3T7|) below. Several other papers (e.g. @ g§ EH 0> S [01 0) 
emphasize, by contrast, the exponential correlation length (= inverse mass gap). The coefficients a\,a2, ■ ■ ■ 
and need not be the same for the two quantities. 



5 This relation between and m/A-^ is based on the result A' 
perturbation theory [^3[ Please note that equation (3) of 
division sign should be multiplication. 



MS /Aj a t = 2 5 / 2 e 7r / 2 ( Ar - 2 ) from lattice 
suffers from a typographical error: the 
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N = 4, 8, however, the agreement with (|1.4j) is better; this suggests that as N increases, 
asymptotic scaling sets in at smaller values of £ [(46] . 

A very similar picture has been obtained by Pade extrapolation of high-temperature series 
PD| , 31, It is found that the critical singularity in the complex /3-plane, which lies at real 
f3 when — 2 < N < 2, bifurcates for N > 2 into a pair of complex-conjugate singularities, 
which tend as N -> oo to = (3/N « 0.32162(1 ± i). In particular, for N = 3, this first 
pair of singularities is located at (3 ~ 1.73 ± 0A2i ||3~2| , Figure 2], extremely close to the real 
axis and to the region in which the Monte Carlo simulations are currently being performed 
((3 = 1.73 corresponds to £ ~ 42). By contrast, for N = 4 the first singularity has moved to 
(3 ~ 2.27 ± O.882, i.e. farther away from the real axis and also towards smaller £ (/3 = 2.27 
corresponds to £ ~ 17). For iV = 8 the first singularity is located at (3 « 3.78 ± 2.38i 
(/3 = 3.78 corresponds to £ ~ 5).[] It is thus very plausible that these complex singularities 
are responsible for preventing asymptotic scaling for £ ~ 100 at N = 3. 

Let us remark that by eliminating (3 from ( |1.4| ) / ( pT5|) , we obtain 



X 



where 



CW£ 2 (log£) 



-(N-l)/(N-2) 



N 



log log £ 
l+cii— t~ + c oi 



log£ 



+ 



dog 2 iog^ 



Cll 



(A^-2) 2 
N-l 



;i.7) 



'1.8a) 



(A^-2) 5 



log 



2tt 



iV-2 



+ (AT - 2) log Q 



+ 



2tt 



N 



■(6i-2oi) (1.8b) 



In finite volume (with periodic boundary conditions), the usual finite-size-scaling Ansatz 
gives, up to corrections of order 



where /e and / x are finite-size-scaling functions, and ^ = 00). Note that lim a 
Q > and likewise for / x . One can also eliminate /3 from ( |1.9|) / ( |1.10|) , yielding 

X (/3,L) = ^,Lf[logai3,L)}^ N -^ N -^f xC (UL) 



(1.9) 
(1.10) 

♦o/e(») = 
(1.11) 



up to corrections of order log log £/ log£. The Ansatze (|1.9|) - (|1.11|) can be written in numer- 
ous equivalent forms, for example, 



X{P,L) 



/3 _ (JV+ l)/(iV-2) e ^/(N~2) f x <£/ L ) 

e(\ogo- {N - imN - 2) m/L) 



1.12) 
;i.i3) 
;i.i4) 



where £ = £(P,L). 



6 We thank Paolo Butera for providing us with precise numerical values of the singularity locations first 
reported in [p2l Figure 2]. 
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2 MGMC Algorithm for Nonlinear cr-Models 



The physical causes of critical slowing-down in the traditional (local) Monte Carlo algo- 
rithms, and the motivation for collective-mode algorithms, have been discussed previously 
0, ||. In particular, paper I |J provides an extensive treatment of the physics behind 
MGMC and the details of its implementation, as well as a comparison with other types of 
collective-mode algorithms (Fourier acceleration and auxiliary- variable algorithms). Here we 
confine ourselves to a brief review of the structure of the MGMC algorithm as applied to a 
nonlinear o"-model. We begin by considering a principal chiral model, i.e. a cr-model taking 
values in some group G C U (N). At the end of this section we shall make some brief remarks 
about the more general case of a cr-model taking values in a homogeneous space M = G/H. 

Consider a principal chiral model with values in G C U(N), defined on a ci-dimensional 
lattice of linear size L; for simplicity let us use a hypercubic lattice with periodic boundary 
conditions. The Hamiltonian is given by 

H{g) = ~P £ Retr(^,) , (2.1) 

(xx 1 ) 

and the Gibbs probability distribution is 

d(i(g) = Z-'expi-Hig)} ]Jdg x , (2.2) 

X 

where dg is Haar measure on G. Our goal is to devise a stochastic updating procedure that 



leaves invariant the probability distribution (2^2) and has an autocorrelation time r which 
is as small as possible. 

The MGMC algorithm for this model is defined as follows |], Sections IV and VII]: 
We first introduce, in addition to the original lattice fl, a sequence of coarse grids Qm = 
fl, f2,/vf-i> . . . , f^o °f linear size L, L/2, L/4, . . . , L/2 M , respectively.^ These coarse grids 

will play a role in intermediate stages of the multi-grid algorithm. Each site y of a grid f^_i 
is considered to be associated with a 2 d block B y in the next-finer grid Qi (Figure |1|). We 
define, therefore, an interpolation operator p\x~\ that maps a field configuration g*-'" 1 ** on 
grid fli-i into one on grid flj, using piecewise- constant interpolation-^ 

(Pi,i-ig {l -% = g { y- 1] for x e B y . (2.3) 

Each coarse-grid field is to be thought of as a correction or update to the next-finer-grid field, 
acting by left multiplication: that is, given a "fine-grid" field g( l < old ) and a "coarse-grid" field 
g^ l ~ l \ we define an updated "fine-grid" field g( l ' new ) = (pi i i^ig( l ~ 1 >)g( l > old ) , or in detail 

g (l,ne W ) = g (l-D g (l,old) for xeBy {2A) 

On each grid Oj, we consider Hamiltonians Hi of the form 

Hfa) = - E Retr(a« , (2.5) 

(xx 1 ) 



7 Clearly L must be a multiple of 2 M . It is therefore most convenient for L to be a power of 2, or at least 
a power of 2 times a small integer. 

8 This update had been proposed earlier, in a "unigrid" context, by Meyer-Ortmanns p%. 
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X 




Figure 1: Staggered coarsening (factor-of-2) in dimension d — 2. Dots (•) are fine-grid 
sites, crosses (x) are coarse-grid sites. 



where a xx , is an iV x iV complex matrix; note that a xx , is not necessarily unitary or a 
multiple of a unitary matrix, and it can vary from one bond (xx') to another. On the finest 
grid (I = M) we have ot^) = (31 for all (xx'} , but on coarser grids more general Hamiltonians 
of the form ( |2.5|) will be generated [see ( |2.7| ) below] . 

For each Hamiltonian Hi of the form (|2.5| ), we assume that there is given a basic stochastic 
iteration g — > Qi(g, Hi) that updates the field configuration g (which lives on grid f^) in such 
a way as to leave invariant the probability distribution ~ exp[—Hi(g)]. For concreteness, we 
shall take Qi to be the single-site heat-bath algorithm for the Hamiltonian Hi. Finally, we 
let 7 be a positive integer, called the cycle control parameter; it will control the number of 
times that the coarse grids are visited. 

The multi-grid Monte Carlo (MGMC) algorithm is then defined recursively as follows: 

procedure mgmc(l, g, Hi) 

comment This algorithm updates a field configuration g in such a 
way as to leave invariant the Gibbs distribution exp[— Hi(g)] dg. 

for j = 1 until m x do g «- Qi(g, Hi) 
if I > then 

compute Hi-t( ■ ) = H t ((pi^i ■ )g) 

for j = 1 until 7 do mgmc(l — 1, h, 

9 *- (Pi,i-ih)g 

endif 

for j = 1 until m<i do g <— Gi(g, Hi) 
end 



We see, then, that a single MGMC iteration at level / consists of three steps: 

1) Pre-sweep. We apply a few (mi) iterations of the basic stochastic iteration (i.e. single- 
site heat-bath) at level I. This produces a new field configuration g that differs significantly 
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from the old one in its short- wavelength components, but not (unfortunately!) in its long- 
wavelength components. 

2) Coarse-grid correction. We compute the coarse-grid Hamiltonian Hi_i(h) = Hi((pii_ih)g). 
As explained earlier, the coarse-grid field h should be thought of as a correction to g. We 
initialize the coarse-grid-correction field h to the identity matrix, and make 7 updates of this 
field using MGMC itself (at level I — 1). We then transfer the updated field h back to grid 

tti using the interpolation operator Pij-i, and use it to correct the field g. The goal of this 
coarse-grid correction step is to update rapidly the long-wavelength components of g. 

3) Post-sweep. We apply, for good measure, a few (1712) more iterations of the basic 
stochastic iteration at level /. 

Let us now look more closely at some aspects of the MGMC algorithm: 

Recursive definition. We have defined the MGMC algorithm recursively — in the coarse- 
grid correction step the procedure mgmc calls itself — because that is by far the simplest 
and clearest approach. However, the MGMC algorithm can be implemented in any com- 
puter programming language, including those (like Fortran) that forbid recursive subroutine 
calls, although the control structure may be a bit more complicated. See [[)8|, Section 4.1] 



for Fortran and Algol implementations of the multi-grid control structure. Our codes are 
available on request (preferably via electronic mail). 



Computation of Assume that the Hamiltonian Hi on level I is of the form (|2.5|) . 

Then a simple computation shows that the coarse-grid Hamiltonian Hi_i(h) = Hi((pi i^\h)g) 

Hi-i(h) = - Retr^^V) + con st , (2.6) 

{yy') 



where 



a y l y' 1] = J2 9x>a£l>gl ■ (2-7) 

(xx') 

x 6 By 

x' 6 B y , 



Note that the coarse-grid Hamiltonian has the same functional form as the "fine- 

grid" Hamiltonian Hf. it is specified by the coefficients {ctyy, 1 ^}- The step ''compute 
therefore means to compute these coefficients. Note also the importance of generalizing the 
Hamiltonian from fl2.1|) to (|2.5| ): even if the fine-grid coefficients are a^} = (31 with (3 > 
independent of (xx 1 ), the coarse-grid coefficients a^ x , for I < M, defined by (|2.7| ), will be 
13 times arbitrary finite sums of elements of G, considered in the space M^v of all complex 
N x N matrices. For G = SU(N) with N > 3, the finite sums of G are dense in M^r. The 
group G = SU(2) is a special case: every linear combination of SU(2) matrices (with real 

9 Note that equation (4.9) of [|| contains a slight error: the sum should be restricted to nearest-neighbor 
pairs (xx') satisfying x S B y , x' £ B y i. A more serious error in || arises from our attempt to allow there 
an arbitrary compact but not-necessarily-unitary group G C GL(N, C). In such a case the nearest-neighbor 
pairs (xx') with x and x' belonging to the same block B y make a non-constant contribution to the coarse- 
grid Hamiltonian, of the form Heti(a y L y ^hyhy). It is thus necessary to generalize the Hamiltonian (|2.5| ) to 
include both nearest-neighbor pairs (xx') and "diagonal" terms x — x'; such a form is then preserved under 
coarsening. We thank Tereza Mendes for pointing out these errors. 
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coefficients) is a nonnegative multiple of an SU(2) matrix, as can be seen from (1.2). Such 



a matrix can therefore be represented in the form ( |1.2j ) but without the constraint <x 2 = 1. 



Note also an important qualitative distinction between the original Hamiltonian (2.1) and 



the "frustrated" Hamiltonian (|2T5|)0: the Hamiltonian ( |2.1| ) has aGxG global symmetry 



g x — > hg x h', while the Hamiltonian (|2.5| ) has only the G symmetry of left multiplication 
g x — > hg x . (If one ng/ii-multiplies g x and g x > by an element h! that does not comute with 
ol^x'i the energy changes.) Fortunately, it is only the left symmetry that is exploited in the 



passage ( |2.4| ) to yet coarser grids. (The coarse-grid-correction move should be a symmetry 
within each block, and contribute an energy only between blocks.) Of course, this fortunate 
circumstance is not a coincidence: from the "unigrid" interpretation |K| || one sees that 
left-multiplication on a very coarse grid acts simply as left-multiplication on a very large 
block on the fine grid — which is of course a symmetry within the block. 

It should also be emphasized that the coefficients in the coarse-grid Hamiltonian (|2.6| ) 
depend implicitly on the current fine-grid field configuration g. Finally, let us emphasize 
that in the MGMC algorithm, unlike the corresponding deterministic MG algorithm, the 
coarse-grid Hamiltonian must (as far as we can tell) be defined by the "variational formula" 
Hi-i(h) = Hi((pij_ih)g)\ only in this way does the algorithm leave invariant the correct 
Gibbs distribution 

Choice of Q\. The basic stochastic iteration Qi( ■ , Hi) can, in principle, be any stochastic 
updating procedure that leaves invariant the probability distribution ~ exp[— Hi(g)]. It 
could, for example, be a single-site (possibly multi-hit) Metropolis algorithm, or a single-site 
quasi-heat-bath |H[ algorithm. Indeed, a quasi-heat-bath algorithm may turn out to be the 



most efficient one in practice (taking into account CPU time). However, for conceptual and 
experimental purposes, a pure heat-bath algorithm is far preferable, as it is a clearly defined 
object with no free parameters. By contrast, the effectiveness of a Metropolis (or quasi- 
heat-bath) algorithm can vary dramatically as a function of the hit size (or underlying Von 
Neumann proposal algorithm) and the Hamiltonian Hi. This is a particularly severe problem 
in MGMC, as the coefficients in Hi are random and highly variable. We have therefore used a 
pure heat-bath algorithm (with red-black ordering) in all our experiments. Technical details 
of this algorithm can be found in Appendix 

Pre-sweeps vs. post-sweeps. The balance between pre-sweeps and post-sweeps is not very 
crucial; only the total nix + m 2 seems to matter much. Indeed, either the pre-sweep or the 
post-sweep (but not both!) could be eliminated entirely. Increasing nix and m 2 improves 
the performance of the MGMC iteration, at the expense of increased computational labor 
per iteration. In this paper we have always taken nix = n^2 = 1, but we make no claim that 
this is optimal. 

Cycle control parameter 7 and computational labor. Clearly, one iteration of the multi- 
grid algorithm at level M comprises one visit to grid M, 7 visits to grid M — 1, 7 2 visits to 



10 The Hamiltonian (2J3) can also be thought of as a principal chiral model in a fixed background "dielectric 
gauge field" 

11 This point has caused some confusion in the literature; for further discussion and references, see |], 
Section X.A]. 
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grid M — 2, and so on. Thus, 7 determines the degree of emphasis placed on the coarse-grid 
updates. (7 = would correspond to the pure heat-bath iteration on the finest grid alone.) 

We can now estimate the computational labor required for one iteration of the multi- 
grid algorithm. Each visit to a given grid involves mi + m 2 heat-bath sweeps on that grid, 
plus the work involved in computing the coarse-grid Hamiltonian and in carrying out the 
interpolation back to the fine-grid. The work involved in all these operations is proportional 
to the number of lattice points on that grid. Let Wi be the work required for these operations 
on grid I. Then, for grids defined by a factor-of-2 coarsening in d dimensions, we have 

Wi ss 2- d{M ~ l) W M , (2.8) 

so that the total work for one multi-grid iteration is 



work(MG) = Y,"1 M ~ lW i 

l=M 

« w M J2(i2~ d ) M ~ l 

l=M 

< W M (l- 7 2- d )- 1 if 7 <2 d . (2.9) 

Thus, provided that 7 < 2 d , the work required for one entire multi-grid iteration is no more 
than (1 — 72~ d )~ 1 times the work required for mi + m 2 heat-bath iterations (plus a little 
auxiliary computation) on the finest grid alone — irrespective of the total number of levels. 
The most common choices are 7 = 1 (the V-cycle) and 7 = 2 (the W-cycle), but we shall 
also consider 7 = 3 (the "super- W-cycle" or S-cycle). 

MGMC vs. renormalization group. We wish to emphasize that although the multi-grid 
method and the block-spin renormalization group (RG) are based on very similar philoso- 
phies — dealing with a single length scale at a time — they are in fact very different. In 
particular, the conditional coarse-grid Hamiltonian H\-\ employed in the MGMC method 
is not the same as the renormalized Hamiltonian given by a block-spin RG transformation. 
The RG transformation computes the marginal, not the conditional, distribution of the block 
means — that is, it integrates over the complementary degrees of freedom, while the MGMC 
method fixes these degrees of freedom at their current (random) values. Our conditional 
Hamiltonian is given by an explicit finite expression, while the marginal (RG) Hamiltonian 
cannot be computed in closed form. The failure to appreciate these distinctions has led to 
some confusion in the literature; for discussion, see || Section X.A]. 

Behavior of Gaussian MGMC. An analogous MGMC algorithm can be devised for the 
Gaussian model (free field) 

H<!P) = ¥>-0 a + §E^> (2.10) 

[xx 1 ) x 

and its behavior can be predicted analytically. A simple Fock-space argument 0, § shows 
that the autocorrelation time of the Gaussian MGMC algorithm is equal to the convergence 
time of the corresponding deterministic MG algorithm for solving Laplace's equation. The 
deterministic MG algorithm is known to behave as follows: 
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Critical slowing-down is completely eliminated for 7 > 1 (i.e. at least a V-cycle) with 
piecewise-linear or smoother interpolation f62fl , and for 7 > 2 (i.e. at least a W-cycle) 
with piecewise- const ant interpolation [p3| . 

For a V-cycle with piecewise-constant interpolation, the exact behavior is not known, 
but it is strongly suspected |j4], section 4.7] and observed in preliminary numerical 



tests [|>f|, |66[ that the critical slowing-down is not eliminated, and that the dynamic 



critical exponent is in fact z — 1. 

The behavior of Gaussian MGMC will be a useful standard of comparison for the behavior 
of MGMC in asymptotically free continuous-spin models. 

Remark. Parisi J74j and Mack and collaborators [ 67], |68 , ^7| |69| have pointed out that 



an update of amplitude A on a block of linear size £b costs an energy AE ~ £ d B ~ 2 A 2 if 
a smooth (i.e. once-continuously-differentiable) kernel is employed, while it costs a much 
larger energy AE ~ l B ~ x A 2 in the case of a piecewise-constant kernel. Thus, the typical 
update amplitude is A ~ £ B in the case of a smooth kernel — and this is the "correct" 

amplitude for an excitation on scale ~ £b in the Gaussian model — while the amplitude 
is only A ~ £ B d in the case of piecewise-constant (step-function) kernel, i.e. a factor 



°b "too small". (This behavior is confirmed numerically fl4"?|.) From these facts one might 



expect that MGMC with piecewise-constant interpolation would not succeed in eliminating 
critical slowing-down for the Gaussian model. However, this conclusion is only partially 
correct: as explained above, the rigorous analysis |63|] shows that piecewise-constant MGMC 



does eliminate critical slowing-down provided that at least a W-cycle is employed (7 > 2). 
Apparently what is happening is that the multiple updates per iteration — £ l B S2j of them 
on length scale £b — compensate for the unduly small amplitude of each individual update. 
Unfortunately, we are unable to make this heuristic argument quantitative, because we do 
not know to what extent these multiple updates are correlated. In particular, we are unable 
to translate the rigorous proof into a heuristic argument in a way that is quantitatively 
correct (i.e. yields the correct critical exponent function of 7). 

Now let us make a few further remarks: 

MGMC for general a-models. Thus far we have considered the MGMC algorithm for 
a o"-model taking values in some group G C U(N). Now let us consider a cx-model taking 
values in some manifold M on which a compact group G C U(N) acts transitively. Then, 
by fixing a reference configuration Tp = {Tp x } x <zQ e M , we can define a (many-to-one) map 
of G n onto M n by 

{9,} — tau • (2.11) 

In particular, Haar measure on G n is mapped onto Haar measure on M Q . Likewise, all 
observables F on M n — including the Hamiltonian — can be "lifted" to G Q : 

F{g) = F(gip) . (2.12) 

Then the expectation value of any observable F(tp) with respect to the Hamiltonian H(ip) 
[and Haar measure on M n ] is equal to the expectation value of F(g) with respect to the 
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Hamiltonian H(g) [and Haar measure on G n }: 

f F{y) e- H ^ I] d V* = I Ho) e~ S(9) U dg x , (2.13) 

where d<p x and dg x are normalized Haar measure on M and G, respectively. Thus, the 
ex-model taking values in M can be "lifted" to a a-model taking values in G. 

Consider, for example, the TV-vector model: here M is the unit sphere in M. N , G is 50(TV) 
or O(N), and the Hamiltonian is 

H (f) = ~ ^2 Ctxx' <fxfx' (2.14) 
(xx') 

where 

&xx' IS cl real number (usually a xx > = (3 for all (xx')). Then the "lifted" Hamiltonian is 

H{g) = - E ^{u xx ,g T x g x ,) , (2.15) 

{xx') 

where 

aw = axx'Wx'W^ (2.16) 

is an x iV real matrix of rank 1. This Hamiltonian is of the form fl2.5|) , and so can be 
handled by the MGMC method described previously. 

In this paper we have chosen to study the 577(2) principal chiral model [= TV- vector model 
with TV = 4] largely because it is easy to implement a heat-bath algorithm on SU(2) ~ S 3 
(see Appendix |A|). The same MGMC principles can be applied to an TV- vector model with 
arbitrary TV, but one then needs to devise a heat-bath algorithm on 50 (TV) or O(TV). 



MGMC with smooth interpolation. Mack and collaborators |67| , |68| , [47], |>£| have advocated 
the use of a smooth interpolation (e.g. piecewise-linear or better) in place of piecewise- 
constant. For fields taking values in a linear state space, this is straightforward to do; but 
for fields taking values in a nonlinear manifold (such as a nonlinear a-model), it is notably 
awkward: how should one define a fractional power of a group element? The solution 



proposed by Hasenbusch, Meyer and Mack [f47| is to define the coarse-grid update by 



g ^ w ) = exp(iA x $>) , (2.17) 

where A x is a suitable smooth kernel and ^ is an element of the Lie algebra of G. (More 
precisely, they take A x to be the longest-wavelength mode of the Laplace operator on an 
x C-b block with Dirichlet boundary conditions.) Note that here the updates are made in 
the "unigrid" style [pOR : the only field in the problem is the original (finest-grid) field; and 
the updates ( |2.17| ) are performed successively on the lxl blocks (in some order), then on 
the 2x2 blocks, then on the 4x4 blocks, etc. For such an update the labor is proportional 



to the number of fine-grid spins involved in the update, i.e. ( |2.8| ) is replaced by 

Wi « W M , (2.18) 
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so that the total work for one unigrid iteration is 

o 

wark{UG) = W M E 7 M "' 

l=M 

MW M if 7 = 1 
j M W M if 7 > 1 

iL d \ogL if 7 = 1 f , 

\L d+l °^ if 7 >l 1 yj 

Thus, the computational complexity per iteration is tolerable only for a V-cycle. Note also 
that the energy H(g( new ^) is an extremely complicated function of the Lie-algebra element 
so it is not feasible to choose \l/ by a heat-bath method; rather, \1/ is chosen by the Metropolis 
prescription with a suitable hit size e ~ £ B . 

It is now natural to ask: what does one gain by the use of a smooth kernel in place of a 
piecewise-constant kernel? In the Gaussian case the answer was given above: with a smooth 
kernel critical slowing-down is eliminated already for a V-cycle[^], while with a piecewise- 
constant kernel a W-cycle is required. However, this fact does not give any advantage to the 
smooth kernel, since in the piecewise-constant case we can readily implement the W-cycle 
in a CPU time of order volume, provided that d > 1 [see (|2.9|) 1. In the non-Gaussian case it 
is less clear what will happen, but we are unable to see why a smooth interpolation should 
handle strongly nonlinear excitations (such as meson-meson scattering or glueball bound 
states) any better than a piecewise-constant interpolation. We return to these points in 
Section O. 



3 Numerical Results 

3.1 Quantities to be Measured 

We begin by defining the quantities — autocorrelation functions and autocorrelation 
times — that characterize the Monte Carlo dynamics. Let A be an observable (i.e. a function 
of the field configuration g). We are interested in the evolution of A in Monte Carlo time, 
and more particularly in the rate at which the system "loses memory" of the past. We define, 
therefore, the unnormalized autocorrelation function^ 

C AA ( t) = (A s A s+t ) - (A) 2 , (3.1) 

12 Here "smooth" refers to an interpolation which is piecewise-linear or better, i.e. one in which linear 
functions are exactly interpolated. It is worth noting that the Hasenbusch-Meyer-Mack kernels |37| |4?| 
do not satisfy this property; indeed, not even a constant function can be exactly interpolated. (That is, the 
constant function cannot be represented as a linear combination of the kernels A x corresponding to non- 
overlapping blocks of the same size, with the exception of the cases £b = 1 and £b = 2.) Therefore, we do not 
know whether critical slowing-down is completely eliminated (for Gaussians) in the Hasenbusch-Meyer-Mack 
V-cycle; we suspect that it may not be. In view of the identities in |], Section VIII], this problem could 
be studied numerically by measuring the convergence rate of the corresponding deterministic MG algorithm 
for solving Laplace's equation. (Note, however, that a more recent work by Hasenbusch and Meyer J6t|] uses 
piecewise-linear interpolation.) 

13 In the mathematics and statistics literature, this is called the autocovariance function. 
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where expectations are taken in equilibrium. The corresponding normalized autocorrelation 
function is 

PA A (t) = C AA (t)/C AA (0) . (3.2) 
We then define the integrated autocorrelation time 



Tint,A — -Z E PAAif) 
A t = -oo 

oo 

= - + E p^(t) 

z t=l 



(3.3) 



[The factor of \ is purely a matter of convention; it is inserted so that T intjA r if Paa(^) ~ 
e -|*l/ r w ith r ^> 1.] The integrated autocorrelation time controls the statistical error in 
Monte Carlo measurements of (A). More precisely, the sample mean 

1 E A (3.4) 



n 



t=i 



has variance 



1 n 

var(A) = — ]T C AA {r-s) 

n r,s=l 

( = — (n— 1) 

1 



(2r ini> A) C AA {0) for n > r 



(3.5) 



(3.6) 



Thus, the variance of A is a factor 2r int ^ A larger than it would be if the {A t } were statistically 
independent. Stated differently, the number of "effectively independent samples" in a run of 
length n is roughly n/2r inttA . The autocorrelation time Tj nt;J 4 (for interesting observables A) 
is therefore a "figure of (de)merit" of a Monte Carlo algorithm. 

We shall measure static quantities (expectations) and dynamic quantities (autocorrelation 
times) for the following observables: 

M 



M' 



E rTx 

X 



(3.7) 
(3.8) 



E 



+ 



Ee 

X 



2nix2/L 



E ^ * ' a x' 

(xx 1 ) 



(3.9) 
(3.10) 



The mean values of these observables give information on different aspects of the 2-point 
function 



G(x) = (cr • cr x ) 
G{p) = E^>o-t 3 



(3.11a) 
(3.11b) 
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In particular, we are interested in the susceptibility 

X = G(0) = V-\M 2 ) 

and the analogous quantity at the smallest nonzero momentum 

F = G(p)\ M=27T/L = V- l {F). 

By combining these we can obtain the (second-moment) correlation length 

\ 1/2 



(MIL 



V4sin 2 (7r/L) J 

(see the Remark below). Finally, we have the (negative) energy 

E = 2G(x)\ M=1 = V-\£) . 
Here V = L 2 is the number of sites in the lattice. 



(3.12) 



(3.13) 



(3.14) 



(3.15) 



Remark. The definition (|3.14j ) is sometimes summarized [42[ by saying that we are 
fitting G(p) to the Ansatz 



G(p) 



r 2 + 4^sin 2 ( Pi /2) 



i=i 



(3.16) 



at p = and \p\ = 2ir/L. This is of course true; but it is important to emphasize that we are 
not assuming that G(p) really has the free-field form ( |3.16|) at general p (of course it doesn't). 
Rather, we simply use this form to motivate one reasonable definition of the second-moment 
correlation length £ on a finite lattice. Another definition — slightly different but equally 
reasonable — would be 



, z(x>/ 

1 x \i=X 



7r) 2 sin 2 (7TXj/L) G(x) 



1/2 



2d 



x) 



4vr 2 



(3.17) 



The two definitions coincide in the infinite- volume limit L — > oo. Finally, let us emphasize 
that neither of these quantities is equal to the exponential correlation length (= 1/mass gap) 



exp 



lim 



— \x\ 



(3.18) 



|x|-»oo \ogG(x) 

However, £ (or £') and £ exp are believed to scale in the same way as (3 — > oo. 

The integrated autocorrelation time r intj A can be estimated by standard procedures of 
statistical time-series analysis fTD), [H]]. These procedures also give statistically valid error 
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d = 2 0(4) model at L = 32 




Type Sweeps Disc 


X S E 


T int,M T int,M 2 


2.00 
2.00 


HB 500000 10000 
W 100000 5000 


92.9 (0.4) 7.61(0.03) 1.15499 (0.00006) 
93.7 (0.3) 7.67(0.02) 1.15525 (0.00010) 


66.19 ( 3.78) 20.68 (0.66) 
1.19 ( 0.02) 2.64 (0.07) 


2.10 
2.10 


HB 500000 10000 
W 100000 5000 


134.6 (0.5) 9.63(0.03) 1.20400 (0.00005) 
134.3 (0.4) 9.61(0.03) 1.20399 (0.00009) 


108.48 ( 7.92) 23.23 (0.79) 
1.33 ( 0.03) 2.74 (0.07) 


2.20 
2.20 


HB 500000 10000 
W 100000 5000 


178.4 (0.6) 11.57(0.04) 1.24789 (0.00005) 

178.5 (0.4) 11.60(0.03) 1.24779 (0.00008) 


145.37 (12.29) 23.97 (0.82) 
1.46 ( 0.03) 2.66 (0.07) 


2.30 
2.30 


HB 500000 10000 
W 100000 5000 


221.2 (0.5) 13.43(0.04) 1.28684 (0.00004) 
221.1 (0.4) 13.41(0.03) 1.28691 (0.00008) 


213.38 (21.88) 20.37 (0.65) 
1.57 ( 0.03) 2.19 (0.05) 


2.40 


W 100000 5000 


258.5 (0.3) 14.94(0.03) 1.32155 (0.00007) 


1.69 ( 0.04) 1.94 (0.04) 



Table 1: Data for the two-dimensional 0(4) model on a 32 x 32 lattice. "Type" is the method 
used (heat-bath, V-cycle, W-cycle or super- W-cycle). "Sweeps" is the number of MGMC 
iterations performed; a heat-bath "iteration" here is equal to two usual heat-bath sweeps. 
"Disc" is the number of iterations discarded prior to beginning the analysis. Standard error 
is shown in parentheses. 

bars on (A) and T intt A- For more details, see JF^, Appendix C]. In this paper we have used 
a self-consistent truncation window of width CT int ^, where c is a constant. In most cases 
we have used c = 6; this choice is reasonable whenever the autocorrelation function pAA(t) 
decays roughly exponentially. However, for the energy (A = £) we noticed that pAA{t) has a 
slowly-decaying component with very small amplitude, and taking c = 6 led to a systematic 
underestimate of Ti n t,s- Therefore, for A = S we took c = 20 for the heat-bath, c = 15 for 
the V-cycle, and c = 10 for the W-cycle and super- W-cycle; this choice seemed to give a 
reasonable window width. 

In setting error bars on £ we have used the triangle inequality; such error bars are overly 
conservative, but we did not feel it was worth the trouble to measure the cross-correlations 
between M 2 and T. 

3.2 Summary of our Runs 

We performed extensive runs using the MGMC algorithm with 7 = 0, 1, 2, 3 (heat-bath, 
V-cycle, W-cycle and S-cycle, respectively) and m\ = m 2 = 1, on lattices of sizes L = 
32, 64, 128, 256. We made an extremely detailed effort to study the W-cycle and a moderately 
detailed effort to study the V-cycle; only a few runs were made to study, for comparison 
purposes, the heat-bath and the S-cycle. In all cases the coarsest grid is taken to be 2 x 2. 
All runs used a random initial configuration ("hot start"). The results of these computations 
are shown in Tables |l], |2|, [| and |. Here (3 = /3o(4) = 2/3su(2) ■ The heat-bath runs use the 
MGMC program but with a cycle control parameter 7 = 0, so each "iteration" consists of 
two heat-bath sweeps. 

In the heat-bath and V-cycle algorithms, the slowest mode (of those we study) is the 
total magnetization Ai: global spin rotations are very slow. In the W-cycle and S-cycle 
algorithms, by contrast, the autocorrelation time r int is negligible: the MGMC algorithm 
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d = 2 0(4) model at L = 64 




Type Sweeps Disc 


X i E 


T int,M T int,M 2 


2.10 
2.10 


V 200000 10000 
W 200000 5000 


161.2 (0.9) 10.45(0.04) 1.20182 (0.00003) 
160.0 (0.6) 10.39(0.03) 1.20187 (0.00003) 


13.47 ( 0.56) 8.49 (0.28) 
1.24 ( 0.02) 3.44 (0.07) 


2.15 


W 200000 5000 


202.9 (0.7) 11.93(0.03) 1.22410 (0.00003) 


1.32 ( 0.02) 3.88 (0.09) 


2.20 
2.20 
2.20 
2.20 


HB 500000 10000 
V 200000 10000 
W 200000 5000 
b 100000 5000 


257.2 (2.3) 13.74(0.09) 1.24532 (0.00002) 

258.2 (1.5) 13.78(0.06) 1.24527 (0.00003) 

258.3 (0.9) 13.79(0.04) 1.24533 (0.00003) 
258.6 (1.1) 13.83(0.04) 1.24534 (0.00004) 


217.47 ( 22.53) 72.96 (4.37) 
18.92 ( 0.93) 11.30 (0.43) 
1.38 ( 0.02) 4.39 (0.10) 
0.51 ( 0.01) 2.89 (0.08) 


2.22 


W 200000 5000 


283.3 (1.0) 14.58(0.04) 1.25349 (0.00003) 


1.42 ( 0.02) 4.58 (0.11) 


2.24 


W 200000 5000 


309.5 (1.0) 15.38(0.04) 1.26149 (0.00003) 


1.45 ( 0.02) 4.58 (0.11) 


2.26 


W 200000 5000 


337.4 (1.1) 16.20(0.04) 1.26932 (0.00003) 


1.48 ( 0.02) 4.71 (0.12) 


2.28 


W 200000 5000 


364.7 (1.2) 16.93(0.04) 1.27692 (0.00003) 


1.51 ( 0.02) 4.79 (0.12) 


2.30 
2.30 
2.30 
2.30 


HB 500000 10000 

\ j c\r\r\r\r\r\ -t aaaa 

V 200000 10000 
W 200000 5000 
b 100000 5000 


393.5 (3.3) 17.71(0.12) 1.28442 (0.00002) 
395.8 (2.0) 17.83(0.07) 1.28441 (0.00003) 
395.3 (1.2) 17.82(0.04) 1.28442 (0.00003) 
394.1 (1.4) 17.79(0.05) 1.28445 (0.00004) 


366.26 ( 49.44) 90.60 (6.04) 
23.84 ( 1.31) 12.39 (0.49) 
1.52 ( 0.02) 5.00 (0.13) 
0.52 ( 0.01) 3.18 (0.09) 


2.35 


W 200000 5000 


471.5 (1.3) 19.87(0.05) 1.30244 (0.00003) 


1.60 ( 0.02) 5.02 (0.13) 


2.40 
2.40 
2.40 
2.40 


HB 500000 10000 
V 200000 10000 
W 200000 5000 
b 100000 5000 


559.1 (3.6) 22.24(0.14) 1.31957 (0.00002) 
549.3 (2.1) 21.86(0.08) 1.31945 (0.00003) 
550.0 (1.3) 21.89(0.05) 1.31942 (0.00002) 
551.9 (1.5) 21.91(0.06) 1.31956 (0.00003) 


545.55 ( 90.06) 93.08 (6.29) 
29.87 ( 1.84) 12.27 (0.49) 
1.66 ( 0.02) 4.83 (0.12) 
0.52 ( 0.01) 3.04 (0.09) 


2.43 


W 200000 5000 


601.5 (1.3) 23.22(0.05) 1.32931 (0.00002) 


1.71 ( 0.03) 4.39 (0.10) 


2.45 


W 200000 5000 


632.5 (1.2) 23.97(0.05) 1.33563 (0.00002) 


1.74 ( 0.03) 4.21 (0.10) 


2.50 
2.50 
2.50 


HB 500000 10000 
V 200000 10000 
W 200000 5000 


706.6 (3.8) 25.79(0.16) 1.35092 (0.00002) 
707.8 (1.9) 25.81(0.08) 1.35092 (0.00002) 
706.4 (1.2) 25.75(0.05) 1.35093 (0.00002) 


713.40 (133.84) 97.46 (6.74) 
35.21 ( 2.29) 10.14 (0.36) 
1.79 ( 0.03) 4.00 (0.09) 


2.55 


W 200000 5000 


779.0 (1.1) 27.49(0.05) 1.36545 (0.00002) 


1.84 ( 0.03) 3.45 (0.07) 


2.60 


W 200000 5000 


848.0 (1.1) 29.08(0.05) 1.37928 (0.00002) 


1.92 ( 0.03) 3.27 (0.07) 



Table 2: Data for the two-dimensional 0(4) model on a 64 x 64 lattice. "Type" is the method 
used (heat-bath, V-cycle, W-cycle or super- W-cycle). "Sweeps" is the number of MGMC 
iterations performed; a heat-bath "iteration" here is equal to two usual heat-bath sweeps. 
"Disc" is the number of iterations discarded prior to beginning the analysis. Standard error 
is shown in parentheses. 
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d = 2 0(4) model at L = 128 


P 


T\rr"\r» ^"YxrAOT^c Tli 


X s a 


T int.M T int,M 2 


2.10 


W 200000 5000 


161.0 ( 0.6) 10.31(0.03) 1.20177 (0.00002) 


1.15 (0.01) 2.55 (0.05) 


O OA 

z.zll 

o on 
z.zU 


a j onnnnn 1 aaaa 
V zUUUUU 1UUUU 

W zUUUUU OUUU 


OCA T / O A\ 1 A 1/?/n AT\ 1 O 4 CAO /A AAAAO\ 

zoy. i ( z.Uj i4.1o(U.U() l.z45Uo (U.UUUUz) 
zoi.i ( i.ij io.yy(^u.U4j i.z4ouy (u.uuuui ) 


oa c a / 1 c\a\ io ca (c\ ca\ 

zU.o4 (1.1)0) Iz.oL) (U.5U) 
l.oU (U.UzJ o.Oi (U.Uoj 


9 9^ 

z.zo 


VV ZUUUUU OUUU 


o4U.o \ l.OJ ID.zD^U.UOJ 1.Z00U4 (U.UUUU1 ) 


l.oJ (U.UZj 4.o4 (U.1UJ 


2.30 
2.30 


V 200000 10000 
W 200000 5000 


452.4 ( 3.9) 19.07(0.11) 1.28394 (0.00001) 
449.0 ( 2.0) 18.83(0.06) 1.28395 (0.00001) 


28.36 (1.70) 18.79 (0.92) 
1.47 (0.02) 5.08 (0.13) 


2.35 


W 200000 5000 


580.0 ( 2.7) 21.81(0.07) 1.30186 (0.00001) 


1.54 (0.02) 6.02 (0.17) 


2.40 
2.40 


V 200000 10000 
W 200000 5000 


747.0 ( 6.8) 25.24(0.16) 1.31885 (0.00001) 

744.1 ( 3.4) 25.12(0.08) 1.31887 (0.00001) 


39.87 (2.84) 26.50 (1.54) 
1.61 (0.02) 6.56 (0.19) 


2.45 


W 200000 5000 


952.9 ( 4.3) 29.08(0.10) 1.33499 (0.00001) 


1.70 (0.03) 7.61 (0.24) 


2.50 
2.50 
2.50 


V 200000 10000 
W 200000 5000 
S 100000 5000 


1172.5 ( 9.8) 32.60(0.21) 1.35028 (0.00001) 
1195.9 ( 4.9) 33.34(0.11) 1.35029 (0.00001) 
1193.5 ( 5.4) 33.21(0.11) 1.35031 (0.00002) 


53.35 (4.40) 30.20 (1.87) 
1.75 (0.03) 7.91 (0.25) 
0.50 (0.01) 4.48 (0.15) 


2.55 


W 200000 5000 


1463.6 ( 5.3) 37.66(0.11) 1.36488 (0.00001) 


1.82 (0.03) 7.78 (0.24) 


2.60 
2.60 


V 200000 10000 
W 200000 5000 


1745.9 (11.1) 41.83(0.24) 1.37874 (0.00001) 
1723.2 ( 5.6) 41.33(0.12) 1.37875 (0.00001) 


68.42 (6.37) 30.01 (1.86) 
1.89 (0.03) 7.72 (0.24) 


2.65 


W 200000 5000 


2003.5 ( 5.6) 45.39(0.12) 1.39195 (0.00001) 


1.95 (0.03) 7.21 (0.22) 


2.70 
2.70 


V 200000 10000 
W 200000 5000 


2290.0 (10.8) 49.53(0.25) 1.40459 (0.00001) 
2294.7 ( 5.2) 49.61(0.12) 1.40460 (0.00001) 


73.51 (7.10) 25.73 (1.47) 
1.99 (0.03) 6.23 (0.18) 


2.75 


W 200000 5000 


2561.3 ( 5.0) 53.26(0.12) 1.41666 (0.00001) 


2.06 (0.03) 5.94 (0.16) 


2.80 


W 200000 5000 


2813.8 ( 4.6) 56.54(0.11) 1.42821 (0.00001) 


2.14 (0.04) 5.12 (0.13) 



Table 3: Data for the two-dimensional 0(4) model on a 128 x 128 lattice. "Type" is 
the method used (heat-bath, V-cycle, W-cycle or super- W-cycle). "Sweeps" is the number 
of MGMC iterations performed; a heat-bath "iteration" here is equal to two usual heat- 
bath sweeps. "Disc" is the number of iterations discarded prior to beginning the analysis. 
Standard error is shown in parentheses. 



18 



d = 2 0(4) model at L = 256 


(3 


Type Sweeps Disc 


X £ E 


T int,M T int,M 2 


2.50 
2.50 


V 200000 10000 
W 100000 5000 


1316.2 (16.6) 34.80(0.30) 1.35021 (0.00001) 
1319.0 (10.2) 35.04(0.19) 1.35021 (0.00001) 


65.79 ( 6.01) 37.72 (2.61) 
1.66 ( 0.04) 7.13 (0.31) 


2.60 
2.60 


V 200000 10000 
W 200000 5000 


2240.7 (30.1) 46.77(0.44) 1.37859 (0.00001) 
2229.4 (12.8) 46.64(0.19) 1.37860 (0.00001) 


86.32 ( 9.04) 53.02 (4.35) 
1.84 ( 0.03) 9.80 (0.34) 


2.70 
2.70 


V 200000 10000 
W 200000 5000 


3648.7 (45.2) 61.76(0.59) 1.40444 (0.00001) 

3646.8 (19.3) 61.92(0.25) 1.40443 (0.00000) 


113.19 (13.59) 62.62 (5.58) 
1.98 ( 0.03) 11.55 (0.44) 


2.80 
2.80 


V 200000 10000 
W 200000 5000 


5519.2 (54.3) 79.39(0.71) 1.42806 (0.00000) 
5448.7 (23.4) 78.33(0.29) 1.42806 (0.00000) 


141.24 (18.91) 64.76 (5.86) 
2.12 ( 0.04) 11.59 (0.44) 


3.00 


W 100000 5000 


9300.0 (30.2) 108.21(0.41) 1.46979 (0.00001) 


2.31 ( 0.06) 8.24 (0.38) 



Table 4: Data for the two-dimensional 0(4) model on a 256 x 256 lattice. "Type" is 
the method used (heat-bath, V-cycle, W-cycle or super- W-cycle). "Sweeps" is the number 
of MGMC iterations performed; a heat-bath "iteration" here is equal to two usual heat- 
bath sweeps. "Disc" is the number of iterations discarded prior to beginning the analysis. 
Standard error is shown in parentheses. 



performs global spin rotations very efficiently. Indeed, if we had taken our coarsest grid to 
be 1 x 1 instead of 2 x 2, then the sequence {M.t\ would be completely uncorrelated, i.e. 

T rat,M = \ fOT a11 

Of greater physical interest is the squared magnetization Ai 2 , which corresponds to 
relative rotations of the spins in different parts of the lattice; this mode is very slow (though 
not the slowest) in the heat-bath and V-cycle algorithms, and it is essentially the slowest 
mode of the W-cycle and S-cycle algorithms. For all these algorithms the autocorrelation 
time r intjM 2 has the same qualitative behavior: as a function of (3 it first rises to a peak and 
then falls; the location of this peak shifts towards (3 = oo as L increases; and the height of 
this peak grows as L increases. In Section |3.4j we shall make a detailed finite-size-scaling 
analysis of these data. 

The autocorrelation time of the energy, r int ) g, is uniformly small: less than about 2.5 for 
the heat bath, and less than about 1 for the other cycles. So we do not bother to report it 
here. 

In Tables ^ and § we summarize, for the convenience of the reader, our best estimates of 
the static quantities \i £ an d E on lattices of size L = 32, 64, 128, 256. These estimates come 
from merging all of the data in Tables |l]-[| together with the data from our Wolff- Swendsen- 
Wang runs |4"2", [EJ-Q These data are consistent with those of Heller |^IJ and Wolff |46[| . We 



use these merged data in our finite-size-scaling analyses whenever x and/or £ is required. 



14 The L = 256, (3 = 2.60, 2.80 data from [Q are unreliable due to an observed metastability, and the 
L = 32, 64, (3 — 2.20, Nun = 1 data from are unreliable due to possibly insufficient equilibration; these 
data points are not included in the means. In all other cases, x, £ and E from the different algorithms agree 
within error bars; this is strong evidence that all the programs are correct! 



19 



Merged 0(4)-Model Static Data 



L 


8 




V 
A 










E 


32 


1.62 


21 


o 


(0 


1) 




("fl 091 

U .U^ 1 


o 


0^41 7 


en 000^41 

lU.UUUu^ 1 


32 


1.70 


27 




(o 


0) 


fi^i 


in m 1 


o 


Q8zt1 S 
yoiio 


(n nnm 1 l 

^U.UUUll J 


32 


1.80 


41 


3 


(0 


0) 


4 fi7 


(n oi l 

U.Ul 1 


1 


nAA^(\ 

U4:4:tJU 


(n nnnm\ 

1U.UUUU 1 1 


32 


1.90 


62 


1 


(0 


1) 


5 Qfi 

O . C*U 


fO 01 1 

U.Ul 1 


1 


1 01 61 

1U1U1 


lo 000071 

1 U. UUUU i I 


32 


2.00 


93 


7 


(0 


1) 


7 68 


fO 01 1 

U.Ul 1 


1 


1 5509 


in 000041 

1 U. UUUU4: 1 


32 


2.10 


134 


G 


(0 


2) 


9 63 


fO 01 1 

U.Ul 1 


1 


90398 


(() 00004 s ! 

1 U. UUUU4: 1 


32 


2.20 


178 


8 


(0 


2) 


1 1 60 

J_ ± . UU 


fO 091 


1 


94787 


(() 000041 

1 U. UUUU4: 1 


32 


2.30 


221 


1 


(0 


3) 


13.42 


In 09) 


1 


9868S 


(() 000041 

1 U. UUUUt: 1 


32 


2.40 


258 


5 


(() 


3) 


14.94 


in 031 

U.UO 1 


1 


OZ;lUlJ 


(() 000071 

1 U. UUUU 1 1 


32 


2.60 


322 


4 


(1 


0) 


1 7 


in i ol 








32 


3.00 


418 


o 


(0 


9) 


21.30 


(0.10) 








32 


3.40 


488 


1 


(0 


8) 


24.30 


(0.10) 








64 


1.70 


27 


9 


(0 


2) 


fi^. 


In 091 








64 


1.80 


41 


2 


(0 


3) 


a fin 


fo i ol 








64 


1.90 


64 


1 


(0 


5) 


r 1 n 


(n i ol 








64 


2.00 


99 


4 


(0 


2) 


7 


(n ni l 

U.Ul J 


1 


luo / o 


(a nnnnzn 

^U.UUUU4: j 


64 


2.10 


160 


5 


(0 


3) 


1 fl 40 

1U.4:U 


(n ni l 

U.Ul 1 


1 


ZiU 1 


1 U.UUUUZi J 


64 


2.15 


202 


9 


(0 


7) 




In o^l 

U.UO J 


1 


99/1 1 n 


^U.UUUUO ) 


64 


2.20 


957 


q 


(() 


4) 


1 3 77 


In 091 

U.UZ J 


I 


9/1 c;qo 
Z4GOZ 


^U.UUUUl ) 


64 


2.22 


983 


■-> 

>J 


( I 


01 

/ 


14.58 


(0.04) 


I 


25349 


(0 00003^ 


64 


2.24 


309 


5 


(1 


0) 


15.38 


(0.04) 


1 


26149 


(0.00003) 


64 


2.26 


337 


4 


(1 


1) 


16.20 


(0.04) 


1 


26932 


(0.00003) 


64 


2.28 


364 


7 


(1 


2) 


16.93 


(0.04) 


1 


27692 


(0.00003) 


64 


2.30 


395 





(o 


7) 


17.81 


(0.03) 


1 


28442 


(0.00001) 


64 


2.35 


471 


5 


(1 


3) 


19.87 


(0.05) 


1 


30244 


(0.00003) 


64 


2.40 


551 


6 


(o 


8) 


21.92 


(0.03) 


1 


31950 


(0.00001) 


64 


2.43 


601 


5 


(1 


3) 


23.22 


(0.05) 


1 


32931 


(0.00002) 


64 


2.45 


632 


5 


(1 


2) 


23.97 


(0.05) 


1 


33563 


(0.00002) 


64 


2.50 


706 


7 


(o 


9) 


25.77 


(0.04) 


1 


35092 


(0.00001) 


64 


2.55 


779 





(1 


1) 


27.49 


(0.05) 


1 


36545 


(0.00002) 


64 


2.60 


847 


9 


(1 


0) 


29.06 


(0.04) 


1 


37928 


(0.00002) 


64 


3.00 


1286 


6 


(3 


o) 


38.40 


(0.20) 








64 


3.40 


1611 


5 


(2 


9) 


45.50 


(0.20) 








64 


3.80 


1865 


7 


(2 


7) 


51.30 


(0.20) 









Table 5: Best estimates of susceptibility, correlation length and energy for the 0(4) model 
on 32 x 32 and 64 x 64 lattices, from Tables |, | and [| 
parentheses. 
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Table 6: Best estimates of susceptibility, correlation length and energy for the 0(4) model 
on 128 x 128 and 256 x 256 lattices, from Tables |, | and (H |9 
in parentheses. 
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3.3 Finite- Size- Scaling Analysis: Static Quantities 



In this subsection and the next, we make our conclusions more quantitative by performing 
a finite-size-scaling analysis. 

We fit L) to the two-loop asymptotic-freedom finite-size-scaling prediction ( |1.12|) by 
plotting £/? 1//2 e _7r/3 versus £/L (Figure |2|). Likewise, we fit L) to the prediction ( [1.1 3D 
by plotting xP 5 ^ e ~ 2nl3 versus £/L (Figure |3|). In both cases the points fall quite nicely 
on a single curve, with the exception of the data points at /3 ^ 2.0 (£ ^8). Note also 
that the scaling curves drop precipitously for £/L > 0.4 (corresponding to ^oo/L > 0.5); 
thus, only those points with £/L < 0.4 bear much relation to the infinite- volume physics. 
By extrapolating the scaling curves to £/L = 0, we can estimate the constants and C x 



defined in O/Q: 

Q(2-ioo P ) = (1.93 ± 0.10) x 10- 2 (3.19a) 
C x(2 -ioop) = (1.88 ± 0.08) x HT 3 (3.19b) 

(subjective 68% confidence intervals). Unfortunately, this value of cannot be directly 
compared with the exact value = 2.349109 . . . x 10~ 2 predicted by the Bethe Ansatz [cf. 
(1.6)1, because we use the second-moment definition of £ [cf. ( |3.14|) 1 while the Bethe- Ansatz 



work |I7|, ^] uses the exponential definition (= inverse mass gap). However, it is found 
empirically that the two definitions of £ agree to within less than 1% — compare our Tables |5] 
and ^| with Table 1 of - - and it would be very interesting to understand why this is so! 

Finally, it should be noted from Tables [| and ^ that the finite-size corrections to £ and 
X at fixed f3 are of order 1-2% when L « 6£ (~ 6^); this finite-size effect should be taken 
into account in future high-precision work. 

The three-loop correction terms a\ and b\ in (1.4) / (|1.5| ) have been computed by Falcioni 
and Treves |^6|: they are 



O! = -—I— [l - (tf - 2)71! - to] (3.20a) 

h = ~ ^ _ 2) [2- 2(N - 2)h - 2h 2 + (N - 1)(1 - k)] (3.20b) 

where Tii, h 2 and k can be computed by evaluating lattice Feynman diagrams; for the nearest- 
neighbor action ( |1 . 1| ) they are h\ ~ —0.09, h 2 ~ 0.52, k ~ 1.57. Thus, for N = 4 we have 
di « —0.053 and b\ ~ 0.031. As already noted by Falcioni and Treves for the case iV = 3, 
these corrections are too small to affect much the quality of the fit in Figures |2| and [| 
However, they do affect the estimates of and C x by a few percent: 

Q(3-ioo P ) ~ C l( 2-ioop) X (l + —z— ) (3.21a) 

Cx(3-ioo P ) ~ C x (2-ioo P ) x ( 1 n — ) (3.21b) 

where (3 av ~ 2.4 is an average value of (3 among those data points contributing to the 
estimates of C% and C x (i.e. L = 128,256 and £/L < 0.2). Hence our best estimates for Q 
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and C x are 

%3-ioop) = (1.97 ± 0.10) x 10- 2 (3.22a) 
C x (3-ioo P ) = (1.86 ± 0.08) x 10" 3 (3.22b) 

(subjective 68% confidence intervals). Unless the higher-order correction coefficients a-i, a^, . . . 
62,^3, • • • turn out to be surprisingly large, the deviations from scaling observed for £ < 8 
will have to be ascribed to nonperturbative effects. 

We also tried to fit x an d £ to the prediction ( |1.14|) by plotting x£~ 2 (log£) 3 / 2 versus 
£/L. However, the fit was rather poor (see Figure f|); and it was not noticeably improved by 
including the universal correction term c\\ log log £/ log £ from ( |1.7|) . Apparently the trouble 
is that | logC^| is comparable in magnitude to log£ for the range of temperatures considered 
here (both are ~ 4), so the correction term coi/log£ is quite significant (and higher-order 
corrections may also be significant). 

In any case, our data show quite good agreement — as good as one has a right to expect 
for this range of /3 — with the two-loop and three-loop asymptotic-freedom predictions. 

Recently, however, Patrascioiu, Seiler and collaborators [153], Q have argued on theo- 
retical grounds that non-Abelian cr-models in two dimensions should have a conventional 
power-law critical point at finite (3 — in sharp contrast to the conventional wisdom ( |L4p - 
( |1.14| ). Moreover, they have claimed to find evidence of such a critical point in a Monte Carlo 
study j55 of the dodecahedron model (a discrete-spin model that is expected to fall in the 



universality class of the 0(4)-symmetric nonlinear cr-model) on lattices of size up to L = 80: 
their best estimates of the critical exponents are v ~ 2.02 and 7 ~ 3.29 (with uncertainty of 
order ±0.1) and thus 2 — rj = j/u ~ 1.627. We therefore attempted to fit our data for £ and 
X to the finite-size-scaling forms 

£{P,L) ~ LfJ((3-/3 c )LV») (3.23) 



X (P,L) ~ f x ({/3 - /3 e )LV v ) (3.24) 
We also fit our data to the Ansatz 

X(P,L) ~ C /V S X ^/L) (3.25) 



[where £ = £(/3,L) and lim^ f x e(x) = C xi > 0], which follows from (ggp/(ggg) by 



eliminating (3. This latter Ansatz does not test the existence of a power-law critical point at 
finite /3, but rather tests the value of the exponent 7/1/ without needing to estimate j3 c . 
First we attempted a fit to (|3.25|) by plotting x£~ 7 ^ versus £/L for a variety of values 



of 7/v. We were unable to obtain a satisfactory fit for any choice of 7/1A In particular, the 
Patrascioiu- Seiler value 7/1/ = 1.627 leads to an extremely poor fit. The least poor fit was 
found at j/u ~ 1.78; it is shown in Figure |5|. However, we do not feel entitled to interpret 
the inadequacy of these fits as evidence against the Ansatz (|3.25|) , because the fit to the 



corresponding asymptotic-freedom Ansatz Ql.14 ) was equally poor 



Next we attempted fits to ( |3.23j ) and ( |3.24j ), which have two and three free parameters, 



respectively. We began by imposing the Patrascioiu- Seiler values v = 2.02 and 7 = 3.29 and 
varying f3 c until we obtained the best fit; the results are shown in Figure ||. We see that at 
f3 < 2.1 — 2.3 (corresponding to £ < 10 — 19), the fits are reasonable. However, there are 
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Figure 4: Finite-size-scaling plot of x£ 2 (log£) 3 / 2 versus £/L, for lattice sizes L = 32 (□), 
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very large deviations at larger values of /3, which cannot be removed simply by varying (3 C . 
Therefore we tried an unbiased study: at each of various trial values of /3 C , we first varied 
v until we obtained the best fit to (|3.23| ); then we fixed this value of v and varied 7 until 
we obtained the best fit to ( p.24| ). Some typical results are shown in Figures [7] ((3 C = 3), [S| 
(P c — 4), I (/3 C = 5), and |10| {(3 C = 10). For the smaller values of j3 c , the fit is still rather 
poor; as /3 C is raised, the fit remains good out to larger values of (3. Moreover, the optimal 
values of v and 7 grow rapidly with fl c : they are in fact very roughly proportional to /3 C (or 
rather, to f3 c — const with const ~ 2). This is exactly the behavior one would expect if the 
asymptotic-freedom Ansatz ( |1.4j )-( |LT4T ) holds, since a power-law critical point at large f3 c 
and with an exponent proportional to (3 C can mimic (for (3 considerably smaller than j3 c ) an 
exponential growth: 




3 C/3 



A + o (I 



(3.26) 



We believe that our data rule out a power-law critical point for /3 C < 5 (a. value well 
beyond the betas at which we ran); and our data are compatible with a power-law critical 
point at (3 C > 5 only if one assumes implausibly large values of v and 7. We conclude that 
the two-dimensional 0(4) model almost certainly does not have a power-law critical point 
at finite /3. 



3.4 Finite-Size-Scaling Analysis: Dynamic Quantities 

Next we apply finite-size scaling to the dynamic quantities r int ^ and r intjM 2. We use 
the Ansatz 

TintAP,L) ~ t{P,Ly^g A {t{(3,L)lL) (3.27) 

for A = Ai,Ai 2 . Here qa is an unknown scaling function, and 5u(0) = linx^o 9a{%) is 
supposed to be finite and nonzero.^ We determine z int ,A by plotting T intjA t;~ Zint ' A versus 
£/L and adjusting z^a until the points fall as closely as possible onto a single curve (with 
priority to the larger L values). We emphasize that the dynamic critical exponent z^a is 
in general different from the exponent z exp associated with the exponential autocorrelation 

Let us begin with r intyM 2. Using the procedure just described, we find 

2.0 ± 0.15 for the heat-bath 

. 1.13 ± 0.11 for the V-cycle , , 

z in t,M 2 ~ S o 60 ± 07 for the W _ cycle i 6 - 2 *) 

0.53 ±0.10 for the S-cycle 



15 It is of course equivalent to use the Ansatz 

UntAPiL) - L Zi ^ A h A (t{p,L)/L) , 

and indeed the two Ansatze are related by hA(x) = x Zint - A g A {x). However, to determine whether 
linizjo g A (x) = lim^^o x~ Zint - A h A (x) is nonzero, it is more convenient to inspect a graph of g A than one 
of h A . 
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Figure 7: Finite-size-scaling plot of £/L (a) and %L 7 /" (b) versus (P-PJL 1 /", for lattice 
sizes L = 32 (□), 64 (x), 128 (O), 256 (O). Here (3 C = 3.0, v = 2.6 and 7 = 4.5. Error bars 
are in most cases invisible. 



30 



•W 1 



0.0 



0(4): ,|8 q = 4.0, ^ = 5.2 




-5 -4 -3 -2 



1 



(a) 



2.0 



0(4): jg =4.0, v = 5.2, 7 = 

"I — M — 1 — I — I — — I — I — I — I — I — I — I — I — I — I — I — 



1.5 



x 



0.5 



0.0 




-5 -4 -3 -2 



1 



(b) 



Figure 8: Finite-size-scaling plot of £/L (a) and 7L 7 /" (b) versus (P-PJL 1 /", for lattice 
sizes L = 32 (□), 64 (x), 128 (O), 256 (O). Here (3 C = 4.0, v = 5.2 and 7 = 8.9. Error bars 
are in most cases invisible. 



31 



1.0 



•W 1 



0.0 



0(4): j8 o = 5.0, ^ = 8.0 




(a) 



1.5 



0(4): B =5.0, iy = 8.0. 7=13.9 

1 - 1 — H — i — i — i — i — i — i — i — i — i — i — i — i — 



1.0 



x 



0.5 



0.0 



(/S-/SJL 1 /" 



(b) 



Figure 9: Finite-size-scaling plot of £/L (a) and xL 7 ^ (b) versus (P-PJL 1 /", for lattice 
sizes L = 32 (□), 64 (x), 128 (O), 256 (O). Here f3 c = 5.0, v = 8.0 and 7 = 13.9. Error 
bars are in most cases invisible. 
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Figure 10: Finite-size-scaling plot of £/L (a) and \L 7 ^ (b) versus (P - p c )L 1/u , for lattice 
sizes L = 32 (□), 64 (x), 128 (O), 256 (O). Here f3 c = 10.0, v = 22.3 and 7 = 39.6. Error 
bars are in most cases invisible. 
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(subjective 68% confidence limits). In Figures |TT]-0 we show the "best" finite-size-scaling 
plots for each case. Note that the corrections to scaling are quite strong: the L = 32 points 
deviate considerably from the asymptotic scaling curve, and even for L = 64 the deviation is 
noticeable (see Figures [T2| and However, it is reasonable to hope that these corrections 
to scaling will have decayed to a small value by L = 128; if so, our estimates of z inttM 2 for 
the V- and W-cycles — which attempt to place the L = 128 and L = 256 points on top of 
one another — will be afflicted by only a small systematic error. In any case, the error bars 



in ( 3.28|) take account of this potential systematic error. 



Note also that the finite-size effects on dynamic quantities are extremely strong, even at 
^/L < 0.2 where the finite-size effects on static quantities are very weak: compare Figures [Tl 
|I"4| with Figures Indeed, in the heat-bath case (Figure [Tip it is far from clear that 



9m 2 (®) — hm^jo 9m 2 ( x ) is finite as it should be; additional data at smaller values of £/L 
would be needed to verify this. In the W-cycle case it does seem likely that Qm 2 {^) is 
finite and nonzero (~ 0.3?); but the scaling function g M 2(x) may well vary by as much as a 
factor of 2 in the range x = £/L < 0.1. We conclude that finite-size corrections to dynamic 
critical behavior can be surprisingly strong; therefore, serious studies of dynamic critical 
phenomena must include a finite-size-scaling analysis. It can be very misleading to assume 
that the finite-size corrections to dynamic quantities are small simply because £/L is small, 
or because the finite-size corrections to static quantities are small. 
Applying the same procedure to r int ^, we obtain 

f 2.0 ± 0.20 for the heat-bath , , 

z int,M ~ \ 1.22 ±0.09 for the V-cycle ^ ' 



(subjective 68% confidence limits). In Figures 15 and 16 we show the "best" finite-size-scaling 



plots for each case. For this observable the finite-size corrections are much weaker. 

3.5 Computational Work 

In Table [7], we show the CPU time per iteration for our MGMC program running on 
one processor of a Cray Y-MP 8/832, as a function of L = 32, 64, 128, 256 and 7 = 0, 1, 2, 3. 
The data deviate from the predicted behavior (|2.9|) for two reasons: Firstly, the fine-grid 



Hamiltonian ( [2.1D is simpler than the coarse-grid Hamiltonians fl2.5|) , and we have exploited 
this by writing special streamlined versions of the heat-bath and "compute i?z_i" subroutines 
for the finest grid alone. Therefore, the CPU time on each grid other than the finest is a 
factor C ~ 2 greater than the estimate ( |2.9| ) assumes. Secondly, the vectorization is more 
effective on the larger lattices because the vectors are longer, so the CPU time grows less 
than proportionately to the volume, at least through L = 256. 
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Figure 11: Finite-size-scaling plot of r inty M 2 i Zint ' M ' 2 versus ^/L for the heat-bath algorithm, 
for lattice sizes L = 32 (□) and 64 (x). Here z intjM 2 = 2.0. 
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Figure 12: Finite-size-scaling plot of r intyM 2^ z ™t,M 2 versus £/L for the V-cycle algorithm, 
for lattice sizes L = 64 (x), 128 (Q) and 256 (O). Here = 1.13. 
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Figure 13: Finite-size-scaling plot of T intyM 2^ Zint ^ 2 versus £/L for the W-cycle algorithm, 
for lattice sizes L = 32 (□), 64 (x), 128 (Q) and 256 (O). Here z m( ^ 2 = 0.60. 
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Figure 14: Finite-size-scaling plot of r intj M 2 C Zint ' M2 versus £/L for the super- W-cycle 
algorithm, for lattice sizes L = 64 (x) and 128 (O)- Here Zi n t t M 2 

= 0.53. 



38 



1.5 



(4) Model; im cycle; z 2 



1.0 



0.5 



□ 



□ >c 



0.0 



0.0 0.1 



0.2 0.3 
f/L 



0.4 0.5 



Figure 15: Finite-size-scaling plot of r int ^ z im,M versus £/L for the heat-bath algorithm, 
for lattice sizes L = 32 (□) and 64 (x). Here z intJ ^ = 2.0. 
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Figure 16: Finite-size-scaling plot of r int ^ z ™t,M ve rsus £/L for the V-cycle algorithm, 
for lattice sizes L = 64 (x), 128 (Q) and 256 (O). Here z int ^ = 1.22. 
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CPU timings for 0(4) MGMC 


L 


7 = 


7 = 1 


7 = 2 


7 = 3 


32 


5.9 


12.2 


33.0 


92.4 


64 


18.7 


33.8 


88.9 


301.2 


128 


66.3 


109.6 


256.7 


986.2 


256 


252.8 


397.6 


815.3 


3305.2 



Table 7: CPU times in milliseconds per iteration for MGMC algorithm for the two- 
dimensional 0(4) model with m\ = m-i = 1 and 7 = 0, 1, 2, 3, on a Cray Y-MP 8/832. 

The running speed for our program at L — 256 was about 100 MFlops. The total CPU 
time for the runs reported here was about 650 hours. 

The statistical efficiency of a Monte Carlo algorithm is inversely proportional to the 
integrated autocorrelation time for the observable (s) of interest, measured in CPU units. 
For the W-cycle at L = 256, (3 = 2.6 (corresponding to £ ~ 47), we have r int ^ M 2 ps 8 Cray 
seconds for a W-cycle, compared to ps 43 seconds for a V-cycle, and ~ 280 seconds for a 
heat bath. (The heat-bath value is an extrapolation from the L = 64, f3 = 2.2 point, using 
finite-size-scaling together with our result z intyM 2 HB ps 2.) Thus, on a 256 x 256 lattice the 
W-cycle MGMC algorithm is about 35 times as efficient as a heat-bath algorithm. 



4 Discussion 

4.1 Heuristic Predictions for z 

In paper I [EL Section IX. C] we predicted that for an asymptotically free theory with a 
continuous symmetry group, such as the O(N) cr-model (N > 2), the MGMC algorithm 
using piecewise-constant interpolation and a W-cycle would completely eliminate critical 
slowing-down except for a possible logarithm. Our data in Section |3]4] show clearly that this 
prediction is incorrect: the exponent zw-cyde = 0.60 ± 0.07 is manifestly not equal to zero. 
It is obviously of great importance to understand heuristically what is going on here. 

Our argument in |6j went as follows: For an asymptotically free theory with a continuous 
symmetry group, the important excitations at short wavelengths are weakly-interacting spin 



waves, i.e. the theory at short wavelengths is almost Gaussian. Now, we know |63[] that 
MGMC with piecewise-constant interpolation and a W-cycle completely eliminates critical 
slowing-down in a Gaussian model. Therefore, we expect that MGMC will almost perfectly 
handle the short-wavelength modes of an asymptotically free theory. On the other hand, 
we know that the long-wavelength modes of an asymptotically free theory are strongly non- 
Gaussian (e.g. there is strong scattering fll4j|), and we do not expect MGMC to handle these 
modes well. Suppose, therefore, that the modes of wavelength < £ max are well handled by 
MGMC moves on the corresponding grids (in the sense that the autocorrelation time for 
these modes is of order 1), but that coarse-grid moves on scales > £ max are ineffective. In 
this case, the correlation length "seen" on the grid of scale £ max will be of order £,/£ ma x', and 
since further coarse grids are ineffective, we expect the autocorrelation time of the MGMC 
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algorithm to be of order tmgmc ~ (£,/£max) ZHB , where z HB is the dynamic critical exponent 
of the single-site heat-bath algorithm. 

The question is: how big is £ max ? In we assumed that £ max ~ £ (for example, 
£max ~ £/10) or at worst £ max ~ £/log p £ - - in which case we would have t mgmc ~ 1 
or tmgmc ~ log 9 £, respectively. But this argument is, in retrospect, much too sloppy. 
An asymptotically free theory is indeed "almost" Gaussian at short length scales, but this 
Gaussianness is approached very slowly: on scale £, the coupling strength (defined e.g. by 
the truncated 4-point function at momentum p ~ l/£) behaves like g(£) ~ 1/ log(£/£). Now, 
suppose that this non- Gaussianness causes MGMC to make small "mistakes''^; and suppose 
that these "mistakes" accumulate (additively) over successive length scales until they reach 
order 1, at which point MGMC becomes ineffective. The length scales (i.e. grids) are indexed 
by log 2 £. This hand-waving argument suggests that £ max should be given by the relation 

/ v^m d{lost) =sC - (41) 

where C is a constant of order 1. It follows that 

£max ~ e , (4.2) 

where 

p = l-e~ c (4.3) 

lies strictly between and 1 (but we are unable to say more than this). Hence we predict 
that zmgmc — V z hb lies strictly between and z HB — a very weak prediction, but one that 
seems at least to be true. 

(We remark that this heuristic argument is very unstable to small perturbations. Suppose, 
for example, that we take instead g(£) ~ l/log p (£/£) with some exponent p ^ 1. Then if 
p < 1 we get £ 

max ~ ^(1)) while if p ^> 1 we get £max ~ So p — 1 is a very delicate 
borderline case. Thus, our heuristic argument, if in fact it is correct, depends crucially on 
the particular measure of "non- Gaussianness" that we have chosen.) 

On the other hand, for the V-cycle the observed exponents z intj M 2 y-cyde = 1-13 ± 0.11 
and z int v _ cyde = 1-22 ± 0.09 are very close to the believed behavior Zv- cy cie = 1 of the 
Gaussian model [see the discussion following (|2.10| )1. This suggests that in V-cycle MGMC 
(with piecewise-constant interpolation) for an asymptotically free model, the slowest modes 
are indeed the Gaussian-like long-wavelength spin waves. It would be interesting to know 
whether zy -cycle for the 0(4) model is exactly equal to 1 (possibly modulo a logarithm), or 
is slightly higher as our data suggest. 



16 By "mistakes" we do not mean to imply a failure to simulate the correct Gibbs measure; of course, by 
construction, the unique invariant measure of the MGMC algorithm is the correct Gibbs measure. Rather, 
we mean a failure to generate independent samples from this Gibbs measure in a time of order 1. 
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4.2 Piecewise-Constant vs. Smooth Interpolation 



As noted at the end of Section |2], Mack and collaborators |67| , p8| (47| , p9| have advocated 
the use of a smooth interpolation (e.g. piecewise-linear or better) in place of piecewise- 
constant. The key question is this: Does the MGMC algorithm have a different dynamic 
critical exponent depending on whether smooth or piecewise-constant interpolation is used? 

In the Gaussian case, the answer is "yes" for a V-cycle but "no" for a W-cycle, as 
discussed in Section |2|. For non-Gaussian asymptotically free models, we conjecture that the 
answer is the same, but to date no full-scale test has been made. Hasenbusch, Meyer and 
Mack |^7J have reported preliminary data for the 0(3) model suggesting that z int ^ M 2 pa 1.2 
for a V-cycle with piecewise-constant interpolation, compared to z inttM 2 pa 0.2 for a V-cycle 
with smooth interpolation. (These authors have not studied a W-cycle, as this would be 
extremely costly in their "unigrid" approach, but it is safe to assume that zw- cyc ie < zy-cyde-) 
For the piecewise-constant V-cycle, this estimate is very close to ours. But the estimate 
for the smooth-interpolation V-cycle (and hence also the smooth-interpolation W-cycle) is 
considerably lower than our estimate z intM 2 pa 0.6 for the piecewise-constant W-cycle - 
suggesting that the smooth interpolation might indeed lead to a smaller dynamic critical 
exponent than piecewise-constant interpolation, even for a W-cycle. However, this estimate 
is based on runs at only four (/3, L) pairs, which moreover have different values of £oc/L, so 
it is very difficult to perform a correct finite-size-scaling analysis. It would be very useful to 
have data at additional (/?, L) pairs, and with higher statistics. 

It should also be noted that since the dynamic critical exponent for MGMC in asymp- 
totically free models is apparently nontrivial (i.e. z^O, 1,2), it is reasonable to expect this 
exponent to be different for each asymptotically free model: for example, in the O(N) mod- 
els the exponent would depend oniV — albeit probably weakly — just as static exponents 
in nontrivial iV-vector models (such as those in d = 3) depend on N. For this reason, it 
would be very useful to have data comparing the two interpolations at the same value of N. 

A very recent preprint by Hasenbusch and Meyer gives preliminary data for the 577(3) 



and CP 3 a-models in two dimensions, using a V-cycle with piecewise-linear interpolation. 
For CP 3 the data are consistent with the dynamic finite-size-scaling Ansatz (|3.27|) with 



z int,M 2 ~ =t 0-1- F° r 577(3), however, the data are very erratic, and do not appear to be 
consistent with fl3.27|) . More data would be very useful here too. 



4.3 Comparison with Other Collective-Mode Algorithms 

It is of great interest to compare the behavior of MGMC with that of Fourier acceleration 
| T\ [75] , [76|, |77J. The two algorithms are very similar in spirit: both are motivated by the 
Gaussian model (for which they eliminate critical slowing-down entirely); both make addi- 
tive updates of a fixed-shape collective mode. The principal difference is the choice of this 
collective mode: for Fourier acceleration it is a sine wave, while for MGMC it is a square 
wave (or triangular wave if one uses piecewise-linear interpolation). Of course, the "correct" 
collective modes in the Gaussian model are sine waves; but the MGMC convergence proof 
|63| shows that square waves are "close enough" (for a W-cycle). There is one additional dif- 



ference: Fourier acceleration is based on small-amplitude updates (i.e. a Langevin or hybrid 
algorithm), while MGMC is based on arbitrary-amplitude updates (a heat-bath algorithm). 
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These close analogies between MGMC and Fourier acceleration suggest that their be- 
havior may be qualitatively similar, in the sense that they work well for the same models 
and work badly for the same models. More precisely, we expect that on length scales where 
the dominant collective modes are very weakly interacting spin waves, both MGMC and 
Fourier acceleration will handle such modes well. On the other hand, on length scales where 
the dominant collective modes are strongly non-Gaussian — either because the spin waves 
are strongly interacting, as in the a-model at t > £ max , or because these modes are in- 
herently discrete, as in the one-component y? 4 model or the high-temperature phase of the 
two-dimensional XY model — we expect that neither MGMC nor Fourier acceleration will 
work particularly well: they will gain either a constant factor in efficiency over the local algo- 
rithms, or perhaps a small reduction in the dynamic critical exponent. (For some nonlinear 
models, Fourier acceleration might work significantly worse than MGMC, if barrier penetra- 
tion by the small-amplitude updates becomes a limiting factor. This difficulty is unrelated 
to critical slowing-down; for example, it occurs already for the double-well (f A model with 
one lattice site.) 

To test this hypothesis, we would like to compare the dynamic critical behavior of MGMC 
and Fourier acceleration for the cx-model. Unfortunately, there is at present very little 
published data on the dynamic critical behavior of Fourier acceleration. Studies of the two- 
dimensional XY model [75|, [76| were discussed in paper II @. The only other published 
test of Fourier acceleration of which we are aware is a study of the two-dimensional SU(3) 



principal chiral model by Dagotto and Kogut |77| . This article reported r int; M 2 an d Tint,e 
for Langevin and hybrid algorithms with and without Fourier acceleration, on lattices of size 
16 2 and 32 2 , for several temperatures in the scaling region. The results indicate qualitatively 
that z < 2 for the hybrid algorithm with Fourier acceleration, but the data are unfortunately 
too crude to allow quantitative conclusions to be drawn. 

The available data are in any case at least roughly consistent with the conjecture that, 
for each model, Fourier acceleration and MGMC have qualitatively similar, and possibly 
even identical, dynamic critical exponents. However, it would be desirable to have more 
precise and comprehensive measurements of the dynamic critical behavior of the Fourier- 
accelerated Langevin and hybrid algorithms, so as to test this conjecture of dynamic (quasi- 
)universality.[] It would also be desirable to have precise measurements of the dynamic 
critical exponent of MGMC for another asymptotically free model, such as the SU(3) prin- 
cipal chiral model, in order to determine whether z varies from one such model to another. 

A very different type of collective-mode algorithm for iV-vector models was proposed 



very recently by Wolff |41]], and independently by Hasenbusch ||43|| . This algorithm is based 
on embedding a field of Ising variables into the iV-vector model, and then simulating the 
induced Ising model by the Swendsen-Wang [80] algorithm (or a variant thereof |4l|). It can 



be argued heuristically |42| that this algorithm is effective at creating long-wavelength spin 
waves, and extensive numerical tests show that critical slowing-down is entirely (or almost 



17 We emphasize that a fair test must use the exact versions of the Langevin |78| and hybrid |7!| algorithms, 
and should tune the step-size in these algorithms (separately for each f3 and L) so as to optimize the 
autocorrelation time. We emphasize also that Fourier acceleration most likely encompasses a continuously 
infinite family of dynamic universality classes, according to the long-distance (low-momentum) properties of 
the acceleration kernel. 
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entirely) eliminated in the two-dimensional iV-vector models for iV = 2, 3, 4 |41 , [14], f|9| .0 
The principal difference between MGMC and Fourier acceleration on the one hand, and 
algorithms of Swendsen-Wang and Wolff type on the other, is how proposals are made 
for updating the long-wavelength modes. MGMC and Fourier acceleration propose additive 
updates of fixed shape (sinusoidal, triangular or square) and variable amplitude, while cluster 
algorithms like Swendsen-Wang and Wolff allow the system, in some sense, to choose its own 
collective modes. 

At present, the Wolff-Swendsen-Wang algorithm is clearly the best algorithm for simu- 
lating N- vector models (at least in an unbroken phase): it is as efficient as MGMC for the 
two-dimensional XY model in the spin-wave regime, it is somewhat more efficient for two- 
dimensional asymptotically free cr-models, and it is far more efficient for the two-dimensional 
XY model in the vortex regime — and it is much simpler to program. On the other hand, 
the MGMC algorithm has a natural extension to cr-models taking values in an arbitrary ho- 
mogeneous space G/H, and at least in principle to lattice gauge theories || Sections IV and 
V]; but it is as yet far from clear whether an efficient Wolff-type embedding algorithm can 



be devised for these models. Indeed, in a separate work |49[ we argue that the generalized 
Wolff-type embedding algorithm can have dynamic critical exponent z C 2 only if the target 
manifold is a sphere (as in the N- vector model), a product of spheres, or the quotient of such 
a space by a discrete group (e.g. real projective space RP 1 ^' 1 ). Finally, Fourier acceleration 
applies to all continuous-spin models, even in the presence of dynamical fermions (which is 
its greatest advantage); but the restriction to small step size is a disadvantage, and the need 
to optimize several tunable parameters makes it more complicated to test and to use. We 
therefore believe that all three classes of collective-mode algorithms — Fourier acceleration, 
MGMC and embedding algorithms — should be developed in parallel, in order to determine 
the conditions under which each one works well or badly. In particular, we advocate that 
a standard test procedure be developed (which autocorrelation times to measure, what sta- 
tistical procedures to use) so that test results obtained by diverse research groups can be 
readily compared. 
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18 For N — 2 this requires that the algorithm also be effective at creating vortex-antivortex pairs; the 
mechanism for this is explained in [Q. 

45 



A Heat-Bath Algorithm for the N- Vector Model 



In implementing the heat-bath step of the MGMC algorithm for the SU(2) principal 
chiral model, we need a subroutine that generates a random matrix U G SU (2) from the 
probability distribution 

dfi(U) = const x exp(tiM*U)dU , (A.l) 

where M is a fixed matrix lying in the linear span of SU (2), and dll is Haar measure. Since 
M = XV with A = (detM) 1 / 2 > and V G SU(2), and U can be parametrized in the form 
(P-. 2|) , it clearly suffices to generate a random unit vector cr = (cr , cr 1 , cr 2 , cr 3 ) G IR 4 from the 
probability distribution 

du(a) = const x exp(/icr°) cicr , (A. 2) 

where cicr is uniform measure on the unit sphere in R 4 (here h = 2A). More generally, we 
can pose this latter problem for unit vectors in R (N > 2); this problem arises also in 
implementing the heat-bath algorithm for an iV-vector model. The problem then divides 
naturally into two parts: 

(a) Generate the component t = cr° with probability density proportional to 

f(t) = e ht (l-t 2 y N - 3 V 2 dt (A.3) 
restricted to the interval [—1, 1]. We can assume without loss of generality that h > 0. 

(b) Generate cr^ = (cr 1 , . . . , a N ~ 1 ) to be a random vector on the sphere of radius (1 — t 2 ) 1 ' 2 
in M N ~ 1 . 

The principal purpose of this Appendix is to discuss three methods for solving problem (a).fj 
At the end we shall make some brief remarks about the much easier problem (b). 

All of our methods are based on the von Neumann rejection technique |]83| , |82|: given a 
function g(t) > f(t), we generate a random variable T with density proportional to g(t) and 
then accept it with probability f(T)/g(T); we keep trying until success. The acceptance 
probability is 

and the expected number of trials is A~ l . Obviously, g must be simple enough that we can 
conveniently generate a random variable with density proportional to g; but subject to this 
constraint we want g to be as similar to / as possible. Our three methods are based on 
different choices of the trial distribution g: 

Method 1 (valid only for iV > 3). g(t) = e ht restricted to [—1, 1]. Such a random 
variable can be generated by transformation: 

T = h- 1 \og[e- h + (e h - e- h )U) , (A.5) 



19 For the special case N = 2, one additional method (rejection from a Gaussian in the angular variables) 
was discussed in the Appendix of JtJ, and several different methods are discussed in [ pl[ . See also |32| pp. 
473-476]. 
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where U is uniform on [0,1]. The proposal is accepted with probability (1 — t 2 )( N 3 ^ 2 . The 
acceptance probability is 



A 



' J V-l ' 
2 , 



N-3 u2 
6N 



h 2 + 0(h 4 



as h 
as h 







(A.6) 
(A.7) 
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2 (^-3)/2 r(^) fc-(^-3)/3 [l +OC/Z- 1 )" 

By construction this algorithm works perfectly for N = 3. For iV > 3 it works decently but 
not perfectly for small h (and deteriorates there when iV is very large), and badly at large 
h. For the case N — 4, this method was proposed previously by Creutz |84| . 

Method 2 (valid only for N > 3). g(t) = e ht [2(l - t)}^ N -' i)/2 restricted to (-oo, 1]. This 
is expressed more conveniently by making the change of variables s — 1 — t. Then the desired 
distribution is 

f( S ) = e -ft- s (AT-3)/2( 2 _ s )(AT-3)/2 (A>g) 

restricted to [0,2], the trial distribution is 

g(s) = e-^(2s) (7V - 3)/2 (A.9) 
on [0, oo), and the proposal is accepted with probability 
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(A.10) 



if < s < 2 

if s > 2 

The trial S is a Gamma^-^-, ran dom variable; for N a positive integer, such a random 
variable can be generated by combining exponential and Gaussian random variables |82| , p. 
405], and this algorithm is efficient if N is not too large. The acceptance probability is 
Method 2 is 

A = (2nh) 1 / 2 e- h I {N ^ )/2 (h) (A.ll) 
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(A.12) 



where I v is a modified Bessel function. Method 2 is extremely bad for small h (especially for 
large N), but it is quite efficient for large h. For the case N — 4, this method was proposed 
previously by several authors |85|, |86], |87 |. 

Method 3. g(t) = e h (l - t 2 )^- 3 )/ 2 restricted to [-1,1]. This is a Beta(^ li ,^ i ) 

distribution, which can be generated from a pair of Gamma^^Y-j variates or in many other 

ways |82], pp. 428-445]. The proposal is accepted with probability e^* -1 ). The acceptance 
probability is 

A = r(^) {h/2y^l 2 e- h I {N ^ )/2 {h) (A.13) 



1 - h + 0(h 2 ) 

2 (^V-3)/2 7r -l/2 r ^|^ h -{N-l)/2 



l + 0{h- 1 ) 



as h 
as h 





oo 



(A.14) 
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This works excellently for small h, and badly for large h. For the case N = 2, this method 
was studied in J7|: the trial variable T is produced by generating a uniform angle 6 G [0, 2n] 
and setting T = cos 9. 

In Table ^] we show the acceptance probability as a function of h for each of the three 
methods, for the case N = 4. (However, the relevant issue in practice is not acceptance 
probability per se, but CPU time; the three methods have slightly different CPU time per 
trial.) Clearly it is necessary to use a combination of methods, in order to obtain a generator 
whose acceptance probability is bounded away from zero uniformly in In our Monte 
Carlo studies of the 0(4) cr-model, we used Method 1 when h < some h cuto ff, and Method 2 
when h > h cuto g. We vectorized the algorithm as follows: Gather all sites with h < h cuto ff, 
and compress them into a single vector; make one trial of Method 1; scatter the "successful" 
T values; gather and recompress the "failures"; and repeat until all sites are successful. 
Then do the same for sites with h > h cuto jj, using Method 2. Empirically we found that for 
runs in the physically interesting range of (3 values (i.e. (3 ~ 2), the CPU time was almost 
completely insensitive to h cuto g in the range 0.5 < h cuto g < 5. We used h cuto g = 1.68474 in 
our production runs. 

Finally, let us make some remarks about the problem of generating a random vector 
(c 1 , . . . , a 1 ^^ 1 ) on the unit sphere of K^ -1 . For N = 3 this is easy: generate a uniform 
random angle 9 G [0, 2tc] and take the cosine and sine. For N = 4 it is also easy: generate 
a uniform random variable a 1 G [—1, 1], and then generate a random vector (a 2 , a 3 ) on the 
circle of radius [1 — (cr 1 ) 2 ] 1 / 2 as just described. (This is the method we employed.) For N > A 
one can generate first a 1 from a Beta(^=^, ^y^) distribution, then a 2 from a Beta(^=^, ^y^) 
distribution scaled by [1 — (cr 1 ) 2 ] 1 / 2 , and so forth. Numerous other methods are described in 



S2, Section V.4]. 



20 In a Monte Carlo application, the h are themselves random variables, and the CPU time is proportional 
to (j4(/i) _1 ) (not (A(h))~ 1 , as we once learned the hard way!). It is thus important to use a generator in 
which A(h) stays away from zero, even for "rare" values of h. This is especially important in the MGMC 
application, in which an extremely wide range of values of h are observed on the different grids. 
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h 


Acceptance Probability 


Method 1 


Method 2 


Method 3 


0.0000 


0.7854 


0.0000 


1.0000 


0.1000 


0.7851 


0.0359 


0.9060 


0.2000 


0.7841 


0.0922 


0.8228 


0.2521 


0.7833 


0.1243 


0.7834 


0.3000 


0.7825 


0.1543 


0.7492 


0.4000 


0.7802 


0.2168 


0.6838 


0.5000 


0.7774 


0.2772 


0.6257 


0.6000 


0.7740 


0.3343 


0.5739 


0.7000 


0.7700 


0.3873 


0.5276 


0.8000 


0.7656 


0.4361 


0.4862 


0.8603 


0.7627 


0.4634 


0.4634 


0.9000 


0.7607 


0.4806 


0.4491 


1.0000 


0.7554 


0.5212 


0.4158 


1.5000 


0.7242 


0.6724 


0.2921 


1.6847 


0.7114 


0.7114 


0.2596 


2.0000 


0.6889 


0.7631 


0.2153 


3.0000 


0.6199 


0.8545 


0.1312 


4.0000 


0.5618 


0.8961 


0.0894 


5.0000 


0.5152 


0.9191 


0.0656 


10.0000 


0.3810 


0.9612 


0.0243 


20.0000 


0.2749 


0.9809 


0.0088 


30.0000 


0.2259 


0.9874 


0.0048 


40.0000 


0.1963 


0.9906 


0.0031 


50.0000 


0.1759 


0.9925 


0.0022 


100.0000 


0.1249 


0.9962 


0.0008 


200.0000 


0.0885 


0.9981 


0.0003 


500.0000 


0.0560 


0.9992 


0.0001 


oo 


0.0000 


1.0000 


0.0000 



Table 8: Acceptance probability for generating a random variable from ( A.3 ) with N = 4, 
using rejection from e ht (method 1), rejection from gamma (method 2) or rejection from beta 
(method 3). 
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